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Abstract 

We consider a delay differential equation that occurs in the study 
of chronic myelogenous leukemia. The equation was investigated nu- 
merically in [8] , [9] . In [6] we studied in detail the stability of the two 
equilibrium points. In the present work, we shortly remind the results 
of [6 and then concentrate on the study of stability of periodic so- 
lutions emerged by Hopf bifurcation from the non-trivial equilibrium 
point. We give the algorithm for approximating a center manifold at 
a typical point of Hopf bifurcation (and an unstable manifold in the 
vicinity of such a point, where such a manifold exists). Then we find 
the normal form of the equation restricted to the center manifold, by 
computing the first Lyapunov coefficient. The normal form allows us 
to establish the stability properties of the periodic solutions occurred 
by Hopf bifurcation. 
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1 The problem 

We consider the delay differential equation that occurs in the study of 
periodic chronic myelogenous leukemia [S], [U] 



i(t) 



1 + x{t) 



x{t) + k /'^^'-'\ . (1) 



Here /3q, n, 5, 7, r are positive parameters and k = 2e~'^'". The unknown 
function, x(-), should also be nonnegative, being a non-dimensional density 
of cells. We do not insist on the significance of the function x{.) or in that 
of the parameters, since these are extensively presented in |9j. Anyway, 
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inasmuch as the equation represents a good model for the periodic chronic 
myelogenous leukemia, both the equilibrium solutions and the periodic so- 
lutions are important, as are their stability properties. 

In the following we use the space B = C([— r, 0],M) (the space of con- 
tinuous, real valued functions defined on [— r, 0] , with the supremum norm, 
denoted by |x|o). Given a function x : [— r, T) i— M, T > and a < t < T, 
we define the function xt £ B hy xt{s) = x{t + s). 

Equation ^ may be written as 

X = h{xt), (2) 

where /i : ;B i— >■ M, and we impose to this equation the initial condition 

xo = (/> G 0. (3) 

In [B] we proved that problem ([2|), ^ has an unique, defined on [— r, oo) 
bounded solution. 



1.1 Equilibrium solutions [S], [U] 

As shown in [5] , the equilibrium points of the problem are 

XI =0, X2 = {^{k-l)-lf'^. 

The second one is acceptable from the biological point of view if and only if 
is strictly positive that is, if and only if 

y(A^-l)-l>0. (4) 

In terms of r, the above inequality may be written as 

1.1/. 5 



r < r^ax ■= In - [1 + 

and since the delay r is positive, the condition (5//3o < 1 follows. 

The biological interpretation of function (3{x) = /3o/(l + x"-) [BJ shows 
that the condition I3{x2) = 6/{k — 1) > should be fulfilled. This is equiv- 
alent to k > 1. 

The linearized equation around one of the equilibrium points is 

z{t) = -[Bi + 6]zit) + kBizit-r), (5) 

with Bi = I3'{x*)x* + /3(x*), x* = xi or x* = X2. 

The eigenvalues of the infinitesimal generator of the semigroup of op- 
erators generated by equation ^ are the solutions of the characteristic 
equation 

\ + 5 + Bi = kBie-^''. (6) 
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2 Stability of equilibrium solution 



Some conclusions on the stability of equilibrium points are presented in 
[8], [9]. The study in [6] has lead us to the following conclusions (obtained 
by using the results in [1]), conclusions that bring some corrections to the 
stability results in [8], [9]. 



2.1 Stability of xi 

In this case Bi = Pq. The equilibrium point xi is stable as long as it 
is the single equilibrium point, that is when ^{k — 1) < 1. When the 
second equilibrium point occurs {^{k — 1) > 1), xi becomes unstable. For 
^(A; — 1) = 1, equation ([6]) admits the solution A = 0. Hence the change of 
stability occurs by traversing the eigenvalue A = 0. 



2.2 Stability of 

We still rely on [6j in this subsection. 

For X2, 

Bi = (3o[n-{n-l)A]/A^ (7) 

where A = f3o{k - l)/5. 

Let us define for future use 

1(1/6 n 
r„ = — ln<^ - - + 1 

We have the following distinct situations. 

I. i?i < that is equivalent to ^{k — 1) > and, in terms of r, 

equivalent to < r < r„. Two subcases are to be considered. 

I. A. Bi < and S + Bi < 0. In this case ReX < for all eigenvalues A 
if and only if |5 + < |A;-Bi| and 

arccos(((^ + ^i)/fcgi) 1 , . 

where ujq is the solution in (0,7r/r) of the equation a;cot(a;r) = —{5 + Bi). 
Alternately, may be defined as 

ooo = lT-H-{S + Bi)r). (9) 

where T : [0,7r) i— > { — oo, 1] is given by 

rp^y) = / ycot(y), y G (0,7r); ^^^^ 

Hence, when the above conditions are satisfied, the equilibrium solution 
X2 is stable. 
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When r|(5 + = 1, the solution X2 is unstable [6]. 
I.B. < and (5 + 5i > 0. 

In this case, Re\ < for all eigenvalues A if and only if 





} (11) 



where ujq is defined as above. When these conditions are satisfied, X2 is 
stable. 

Remarks, a) If Si < 0, 6 + Bi = 0, the solution is stable if and only 
if —kBir < Tr/2 while, for this case, the point kBir = — 7r/2 is a Hopf 
bifurcation point. 

b) Assume that we vary r and keep all the other parameters fixed. The 
conditions ([5]) or for r are not as simple as they seem, because Bi is 
itself a function of r, (being a function of k). Let us consider the function 



If for a certain r* we have g{r*) = (that is the condition for the change 
of stability), in order to find whether a value ri in a neighborhood of r* is 
in the stability zone or not, we have to know the sign of g{ri), hence we 
have to study the monotony properties of function g' in a neighborhood of 
r*. This will be done in the last section for an example. 

II. Bi > that is equivalent to ^(A: — 1) < and, in terms of r, 
equivalent to r„ < r < r^ax- 

In this case, we can only have 6 + Bi > 0, and, by [T], ReX < for all 
eigenvalues A if and only if kBi < 6 + Bi. 

By using ([7]) we see the above inequality is equivalent to ^"^^"'''^ > 1 that 
is already satisfied, since X2 exists and is positive. It follows that if i?i > 
then X2 is stable. 

3 Hopf bifurcation points 

By denoting A = /x + iw, and by equating the real, respectively the 
imaginary part of the characteristic equation ([6|), we obtain 



g{r)=T ^(— ((5 + i?i (r))r) — arccos 




H + 5 + Bi = fc^ie"'"' cos(a;r) 



(13) 



UJ = —kBie ^'"sin(cjr). 
It is useful to consider the case /u = in the above equations. 



5 + Bi = kBi cos{ujr) 



(14) 



Lo = —kBi sin(a;r). 
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In the above stability discussion, in I., the case 

wccos{{5 + Bi)/{kBi)) 
r = (15j 

Wo 

occurs on the frontier of the stabihty domain. This relation, together with 
a;ocot(wor) = —{6 + Bi) imply that ujq = {kB\)^ — (<5 + B\)^ and that 
the pair //* = 0, w* = wo represents a solution of (fT3]) . 

Hence a point in the parameter space, where relation (|15p is satisfied, 
might be a Hopf bifurcation point when one of the parameters is varied. The 
Hopf bifurcation may be either a non-degenerated or a degenerated one. 

To be more precise, we consider the vector of parameters (/3o, n, 5, 7, r) 
and denote it by a = (aj)i<j<5. 

We choose a point in the parameter space, a* = (a*)i<j<5, such that, 
for this choice of parameters, there are two pure imaginary eigenvalues 
'^i,2(a*) = ±^w*, w* := w(a*) > and all other eigenvalues A have strictly 
negative real parts (i.e. a point where condition (|15p holds). 

In the following we keep fixed all parameters excepting one, aj, that we 
vary (j fixed, 1 < j < 5,). For further simplicity, we denote aj := a, such 
that Ai^2(a) becomes Ai^2(a) and such that ReXi^2{ct*) = 0. 

If ^g^-^' (a*) / 0, i = 1,2, then at a* a Hopf bifurcation occurs and, 
in a neighborhood U* of a*, either to the left or to the right of a* we 
have -ReAi_2(a) < 0, while i2eAi_2(a) > at the opposite side of a*. The 
neighborhood U* is such that for all a £ U* all other eigenvalues besides 
Ai^2 have real parts less than a fixed negative value. 

4 On the stability properties of the periodic solu- 
tions emerged by Hopf bifurcation 

We consider a periodic solution occured by Hopf bifurcation, as described 
in the previous section. In order to establish whether or not such a periodic 
solution is stable, we have to compute the normal form of the reduced to the 
center manifold (or to the unstable manifold) problem. In this section we 
present the algorithm for obtaining this normal form and the consequences 
on the stability of periodic solutions. 

Since we study local bifurcations around X2, throughout this section we 
work with eq. ([1]) translated by X2, that is, the new unknown function 
z = X — X2 is considered. 

4.1 The delay equation as an equation in a Banach space 

We assume that we are in one of the situations, mentioned in Section 3, 
where a Hopf bifurcation occurs. In order to construct the center (or the 
unstable) manifold at the equilibrium point X2 we need to write our equation 
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as a differential equation in a Banach space. We use the method of Teresa 
Faria, from [2], as we did also in [3], [5]. We consider the space 

i3o = s ^ : [— r, 0] iH> M, is continuous on [— r, 0) A 3 lim ^(s) S M > , 

that is formed of functions tp = + ado, where (/9 € a € M and 
do : [-r, 0] ^ M, 

. = 0. 

The space Bo is normed by H-^H = |(/9|o + |cr|. The linearized equation ([5]) can 
be written as 

z = Lzt, (16) 
= 1° dviOMO), 
(the integral is a Stieljes one) with 

{-kBi, s = -r; 
0, sG(-r,0); 
-(Si +(5), s = 0. 

Nonlinear equation ([1]) may be written as 

i{t) = L{zt) + f{z{t),z{t-r)), (17) 

where /(■, •) comprises the nonlinear terms of the equation. Since z{t) and 
z{t — r) are functions of Zf , eq. pT|) may be written as 

z{t) = L{zt) + f{zt). (18) 

In [2], the linear operator A : Ci([-r, 0],M) C Bo ^ Bo 

A{ip) = ^ + do[L{ip) - (19) 

is defined and it is proved that this is the infinitesimal generator of the 
semigroup of operators {S{t)}t>o given by S{t){(p) = zt{(j)), where z{t,4') is 
the solution of equation ()16p with the initial condition zq = Then the 
nonlinear equation may be written as an equation in Bo, that is 

^ = Aizt) + dof{zt). (20) 
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4.2 Space of eigenfunctions corresponding to A12, the pro- 
jector on this space 

We return for the moment to the hnear equation ([5]). The eigenfunc- 
tions corresponding to the two eigenvalues Ai^2 = /"(a) ± iu^a) are given by 
^i,2{s) = e^^'^'', s G [— 0]. Since the eigenfunctions are complex functions, 
we need to use the complexification of the spaces B, Bq that we denote by 
Be, Bqc, respectively. We denote by ^A the subspace of Be generated by 

In [2], for an equation of the form (j20p . a projector defined on Ai is 
constructed. For its construction we need a bilinear form that is defined by 
using the transposed equation of the linearized equation. The transposed 
equation is 

y{s) = - [ yis-9)d[rj{9)], 



hence 

y{s) = {Bi + 5)y{s) - kBiy{s + r). 
The characteristic equation for the adjoint problem is 

A = Bi+5-kBie^\ 

that admits the solutions A1.2 = — /^(a) i ioj{a). The corresponding eigen- 
functions are Vi(C) = e(-^+*'^)^, ^2(0 = e^-^-*^)'^, C G [0,r]. 

In order to construct the projector defined on Bqc with values on M., 
we define the bilinear form [3] 

(V, ^) = mm - f j' i^ic - o)dmmdc = 



= mm+kBi J i^iC + rMOdC, 

where ■0 '■ [0, r] — )■ M. 

We look for linear combinations of the vectors ipj, that we will denote 
^i, i = 1, 2, such that {^t, (pj) = 5ij. For this we determine the 2x2 matrix 
E, with elements Cij = {tpi, (pj). 

en = = MO)MO) + kBi / g^-^+^'^^^'^+^^e^'^+^'^^^d^ = 1 - e'^^ 

J —r 







621 = (^2, ^x) = V2(0)v?i(0)+A;Si / e^-i^~i^){Q^r)^(^.^^^)e^Q ^ x^{b^B^^iiMuj)r, 
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622 = (^-2, ^2) = V'2(0)V^2(0) + kB^ / e^-f^-'^W+r) ^{^.-iu)e ^ ^ _ ^/.r^ 



We have det^ = 611622-612621 = {I - e^''' f -[{1 + {5 + Bi + ii)rf + uj'^r\ 
Since 



, ^2 7 V V'2 

we obtain 

1 — C\2 

*1 = , ^ [622^1 - ei2'(/'2] = -—-^'4)2, ^'2 = ^1- 

det E det E 

We can now write the projector defined on Bqc £iiid with values in Ai. 
It is given, for = + d^a E ;Boc by 

P(V) = ((^1, </>) + ^i(O)a) (^1 + ((^2, 4>) + ^2(0)a) (^2. (21) 

If 0" = and, thus, tp = cj) G Be, we have P(0) = (^'i,(/>)v3i + (^2i'/')<^2- 
Now, for the solution z (-,</>) of (f20|) we can write 

= (/'lUl(t) + 02^2(0 + v(t), 
with ni(i) = (^1, Zt), U2{t) = {^2, Zt), V = (/ - V)zt. 

We project now the equation (PU|) by to obtain, 

^ = 5u + ^(O)7(0iUi+</.2n2 + v), (22) 



where 



^^2 y ' \ A2 / ' V *2 



Remark that ui,U2 are complex functions, but, since our initial problem 
is formulated in terms of real functions, we have ui = U2, and the two scalar 
equations comprised in ()22p are complex conjugated one to the other. Hence 
it suffices to study one of the them, let us say the equation for ni. We denote 
ui = u and the equation to be studied in the following is 

du — 

— = {fi + iuj)u + *i(0)/(^iu + if2U + v). (23) 
at 

For a = a* we have /^(ct*) = hence 

^'^^^ - [l + (5 + i?i)rP+a;V2- ^^^^ 
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4.3 The equation restricted to the invariant (center or un- 
stable) manifold 

For those a €z U* where ReXi^2{ce) ^ 0, since ah other eigenvalues have 
negative real part, there is a local invariant manifold, the local unstable 
manifold, Wf^^, for ReXi^2{ct) > 0, or the local center manifold, Wf^^, at a* 
(it:eAi,2(a*) = 0). 

The local invariant manifold is an at least invariant manifold, tangent 
to the space A4 at the point z = (that is x = X2), and it is the graph 
of a function w{a){-) defined on a neighborhood of zero in Ai and taking 
values in (/ — V)Bc- A point on the local invariant manifold has the form 
u ifi + u ip2 + w{a){u +U(p2)- Since it is an invariant manifold, if the initial 
condition is taken on the manifold, then its image through the semigroup 
{f (t)}t>o, f{t){(t)) = zt{4>), is still on the manifold (?(*)(</)) = T{t){4>) - X2, 
where T{t) is defined at the end of Subsection 1.1). Hence 

f{t){^) = u{t)ipi + u{i)^2 + w{a){u{t)ipi + Mt)v2), (25) 

with n(-), solution of the equation 

du ^ 

— = Ai(a)u + 'I'i(0)/(v3iu + (/52^i + i«(a)(nv3i +M9J2)), (26) 

with the initial condition n(0) = uq, where Vi^i^^ = UQip\^UQi~p2- The real and 
the imaginary parts of this complex equation, represent the two-dimensional 
restricted to the center manifold problem. We can study this problem with 
the tools of planar dynamical systems theory (see, e.g. [7]). 

4.4 The normal form of the reduced equation 

In order to find the normal form of equation (|26|) we must consider the 
series of powers of li, u for w(a) and for f{zt). We set 

w{a){u,u) = ^ -T^Wjk{a)u^u'', (27) 

where ■w{a){u,u) := 'w{a){u{t)ipi +u{t)if2), Wjk{a) G B, u: [— r, 00) i-)- C. 

f{zt) is, due to ([25]) . a function of u, u, and we develop it in a series of 
powers 

^(^*)= E -^/^-fcH"'^'- (28) 
Equation ()26p becomes 

^ = Ai(a)n+ -l-g.^(a)u^u\ (29) 
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where 

gjk{a) = ^'i(0)/,fc(a). (30) 
If the first Lyapunov coefficient at a*, [7], 

^i(a*) = i^^Re{ig2Q{a*)gii{a*) +uj*g2i{a*)) (31) 

is not equal to zero, then at a = a* our equation presents a non-degenerated 
Hopf bifurcation. 

We assume that the initial function, (/>, is on the invariant manifold. We 
give the algorithm for computing the first Lyapunov coefficient. We are 
interested in the nonlinear part of our equation. We consider the function 
(3{x)x, compute its derivatives of order higher than one at the point X2, and 
denote them by 

5n = /3^")(x2)a:2 + n/3("-i)(x2). 

Let F be the function F{z) = ^B2Z^ + ^B^z^ + ■^B^z'^ + .... The nonlinear 
part of the RHS of (HID is J{z{t), z{t - r)) = -F{z{t)) + kF{z{t - r)), or, in 
terms of zt, the function / of ([l8]) is f{zt) = -F{zt{0)) + kF{zt{-r)). 

The normal form of the restriction of the problem to the unstable man- 
ifold is [7] 

— u = (/3(a) + i)n -|- su|np, (32) 

where /3(a) = /u(a)/w(a), s = sign{li{a*)). For a = a*, /3(a*) = 0, and the 
normal form on the center manifold is obtained. The existence and stability 
properties of the cycle emerged by Hopf bifurcation are given by the signs 
of /3(a) (for a G C/*) and /i(a*). We must compute /i(a*), given by (f3T]l . 

In order to identify the coefficients gjk{ot*) we must find the coefficients 
for the series of /. The following computations are done at a = a*, but for 
simplicity, the parameter will not be written. 

We have (by denoting tpi = Lp) 

j{zt) = -F{[fmm + kF{[fm{-r)) = (33) 

= -^S2[uv3(0) + u^(0) + ^u;2o(0)u2 + u;ii(0)nn + ^t«02(0)u^ + ...]^- 
-^S3[uy5(0) + n^(0) + ]^W2q{Q)u^ + wii{{))uu + ^U'02(0)m2 + ...]3 - ....+ 
+ ^A;i?2 [u^{-r) + U(p{-r) + ^W2o{-r)u^ + wii{-r)uu+ ^Wo2{-r)lf + + 
+-^kB3[u(p{-r)+u^{-r)+^W2oi-r)u^+wu{-r)uu+^wo2{-r^^ = 



jlkl-"' 

j+k>2 •' 
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The second order terms are 

-^B2{u^ + 2uu+lf) + + 2™ + nV(-0), 

and the coefficients 

/20 = B2{k^\-r) - 1) = -B2{1 - ke~^^*n, 

hi = B2{k - 1), 

/02 = B2{k^\-r) - 1) = -B2{1 - A:e2-*-), 

follow. Prom these relations and ([30|) . we already find 520 > du, S02- 

In order to determine li we still need to compute the coefficient 521- 
From (1331) we obtain 



/21 = B2 (-2u;ii(0) - u;2o(0) + 2ke-''^'wu{-r) + ke''^'' W2Q{-r)) - B^^ (l - ke 

(34) 

hence we have to determine W2o(0), W2o(— ?'), ifii(O), ?i;ii(— r). 

The functions G are determined from the differential equations 
obtained by equating the same degree terms of ( [TU] , [Ij ) : 

^ = 7]^*fe^'^Vi(5)+ (35) 

j+k>2-'' ' j+k>2-'' 

+ Yl 7^^j-^'^'''''^2(s) + ^ ^ ^u;,fc(s)^/%^ 

j+k>2-'' ' j+k>2-'' 

and the integration constants are determined from0 

E E -n:[5,fc^'^'v'i(o)+ -i^5,,u%V2(o) = 



j+fc>2 i+A:>2 j+k>2 •' 



(36) 



j+k>2-^' ' j+k>2-'' ' j+fc>2 

Hence, the function W2o(") is solution of the differential equation, 
^o(s) = 2^*^20(5) +g2w{s) +9o2^{s)- 



^In the first form of this article, published in Journal of Middle Volga Mathematical 
Society, (journal formerly named Proceedings of ...), unhappily, relation (|36|) was written 
in an erroneous form. The subsequent computations were affected by the error therein. 
The small corrections that we brought to our work in its second version posted on arXiv 
did not cover this error. We make here the necessary corrections, and we apologize to 
the readers that might have been confused by our previous mistakes. An Erratum will be 
send to Journal of Middle Volga Math. Soc. 
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We integrate between —r and and obtain 



W2o{0) - W2o{-r)e^ 



By equating the coefficients of in ([36|) . the relation 

2uj*iw2o{0) +520 +9o2 = -{Bi + 5)w2o(0) + kBiW2o{-r) + /20 
results. From the two relations above, we get 

""^20(0) = c 



e J20 + 920 I — I — e — e | + 



to" 



+902 



'W2o{-r) = c 



kBii ^ kBii ^^^^*^ ^ ^2iuj*r 



3uj* Soj* 
/20 + ff2o(l-2e^ 



-^^02(l + 2e3-*- + ^i(l-e3-*-)) 



where c = [-kBi + {Bi + (5) cos(2w*r) - 2a;* sin(2a;*r) - i(2a;* cos(2a;*r) + 
{Bi + 6) sm{2uj*r))]/[{kBi)^ + {Bi + 6f + (2w*)2 - 2fe5i(5i + (5) cos(2a;*r) + 
AkBiLO* sm{2uj*r)]. 

We now compute wu. It is the solution of 

By integrating between —r and 0, we obtain 



wn{0) = wn{-r) - —511(1 - e"" + —511 (1 - 
u; 

The second relation, is here 

511+911 = -(^1 + S)wn{0) + kBiWn{-r) + fn- 
It follows that 

kBii 



wu{0) = ci 



wn{-r) = ci 



/ii - 511 - ffii + — ^ 911(1 - e-'^ - 511 (1 - e'^ 



fii - 911 - 9u + ^^1±^ ( 511(1 - e-'-) - 511 (1 - e-'-) 



with ci = 1/(^1 + S — kBi). Now we have all the elements for determining 
/21 by formula ((Ml), and then we can compute 



521 = ^'l(0)/21, 
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and after that, the first Lyapunov coefficient, that determines the stability 
properties of the periodic orbit emerged by Hopf bifurcation. 

The above indicated computations, are easy to perform with a computing 
program (Maple e.g.), for a given Hopf bifurcation point, i.e. if all the 
parameters are known. If we try to make the computations with all (or 
some of) the parameters non-determined, we are lead to very complicated 
expressions that can not be handled even by the program Maple itself. 

4.5 An example 

In order to provide an example, we choose the following "strategy" for 
finding Hopf bifurcation points. We choose n*, /3q, 5* ,k* , in an acceptable 
from biological point of view zone of parameters. With these values of 
the parameters we can compute B^, (the value of Bi at these parameters 
values), p* = 6* + B*, q* = k*Bl, and determine u* and r* such that a 
Hopf bifurcation actually takes place. That is, we set w* = \/ {q*)"^ — (p*)^ 
and r* = aiccos{p* /q*)/uj* . Once r* determined, we can also compute 7*. 

As before, we choose one of the parameters aj to be the variable one and 
check the condition that the curve Xi{aj) actually crosses the imaginary 
axis in the complex plane, when Uj varies around a^. This can be done by 
example by differentiating relations ()13p with respect to aj and by solving 
the obtained system of equations with respect to /i'(ap, uj'{a*). 

The example. We consider the following selection of parameters: 

n* = 12, = 1.77, 6* = 0.05, k* = 1.180746972 7V* = 0.527). 

These parameters were chosen in the zone of biological interest, as can be 
seen in [S]. 

We find then X2 = 1.150859618, B^ = -2.524 121872, p* = 6* + B^ = 
-2.474121872, q* = k*Bl = -2.980349005, w* = ^/q*"^ - p*^ = 1.661686238, 
r* = arccosip* /q*)/uj* = .3559207407, 7* = 1.48067. 

We choose to vary r and, by differentiating the relations (jl3p with respect 
to r find that fi'{r*) = 25.66. 

The first Lyapunov coefficient is easily computed since all the parameters 
are known, by following the algorithm in the previous subsection. We get 

li(r*) = -43.71063. 

Thus the normal form of the restriction of the equation to the unstable and 
to the center manifold is (by P2|) ) 

—u = (/3(r) + i)u — n|up, 

with /3(r) = and /3(r*) = 0. The normal form in polar coordinates is 

r p=(/3(r)-p2)p, 
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We are interested in the behavior of the solutions on a small neighborhood 
of r*. Since fi'(r*) > 0, on such a neighborhood r < r* ^ fj,{r) < ^ 
/3(r) < 0, and the equation in p above has the only equilibrium solution 
p = 4^ X = X2 and it is stable, while r > r* =^> p{r) > /3(r) > 0. 
Thus a stable periodic solution exists for r > r* . 




a) r=0.35. 




b) r=0.35592. 




c) r=0.36. 

Fig.l) Trajectories obtained by numerical integration of eq. ([T]) for initial 
condition cl){s) = cos(^s), for some values of r around the Hopf bifurcation 

value r* = .3559207407. 

These assertions seem in contradiction with the condition of stability dH). 
Actually there is no contradiction. The function g defined in (|12|) satisfies, 
for the chosen parameters values, g{r*) = and ^{r*) = —20.236. Thus, 
r < r* ^ g{r) > and the equilibrium point X2 is stable, while for r > r*, 
g{r) < and the equilibrium point is unstable. 

The numerical integration confirms our theoretical analysis, as the figures 
above show. 
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